# usage: just_islands.pl bam mindepth 

unless(-r $ARGV[0]) {
    die "no bam file you loser\n";
}
unless($ARGV[1]) {
    die "no mindepth you loser\n";
}

@islands = ();
$last_start = -1;
$last_pos = -1;
$last_chr = "NULL";
$ok = 0;
(open(DEPTH, "samtools depth $ARGV[0] 2> /dev/null |")) || die "no open\n";
while (<DEPTH>) {
    chomp;
    @fields = split ("\t", $_);

    if($fields[1] != (1 + $last_pos)) {
	if ($ok) {
	    $island = "$last_chr" . "\t" . "$last_start" . "\t" . "$last_pos";
	    push(@islands,$island);
	    $ok = 0;
	}
	$last_start = $fields[1];
    }
    if($fields[2] >= $ARGV[1]) {
	$ok = 1;
    }
    $last_pos = $fields[1];
    $last_chr = $fields[0];
}
if($ok) {
    $island = "$last_chr" . "\t" . "$last_start" . "\t" . "$last_pos";
    push(@islands,$island);
}

foreach $island (@islands) {
    # print what you have so far
    print "$island\t";
    # print the locus size and get the sam query
    @fields = split ("\t", $island);
    $locus_size = $fields[2] - $fields[1] + 1;
    print "$locus_size\t";
    $query = "$fields[0]" . ":" . "$fields[1]" . "-" . "$fields[2]";
    # parse the reads
    %read_sizes = ();
    $total = 0;
    $ok_total = 0;
    (open(SAM, "samtools view $ARGV[0] $query |")) || die "no samtools view open for $ARGV[0] $query\n";
    while (<SAM>) {
	@sf = split ("\t", $_);
	$size = length $sf[9];
	++$read_sizes{$size};
	++$total;
	if(($size >= 20) and ($size <= 24)) {
	    ++$ok_total;
	}
    }
    close SAM;
    @sorted = sort { $read_sizes{$b} <=> $read_sizes{$a} } keys %read_sizes;
    if(($ok_total / $total) >= 0.8) {
	$call = $sorted[0];
    } else {
	$call = "N";
    }
    print "$total\t$call\n";
}

